## ----setup, include=FALSE-------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(data.table)
library(fixest)
library(modelsummary)
library(ggthemr)
options(scipen = 999)
setFixest_nthreads(14)
setDTthreads(14)


## ----Params, echo = FALSE-------------------------------------------------------
FigFolder <- ""


## ----figure_functions-----------------------------------------------------------
coef_to_plot <- function(df, label){
  DT <- as.data.frame(df)%>%
    rownames_to_column(var = "Var") %>% 
    rename( se = `Std. Error`) %>%
    mutate( Var = str_remove(Var, ":MaxTim")) %>%
    mutate( var_n = as.numeric(str_remove(Var, "time_to_treatment::")),
            model = label,
            var_n = replace(var_n, var_n==-13, -1),
            se    = replace(se, var_n==-1, 0),
            Estimate   = replace(Estimate, var_n==-1, 0),
            ) %>%
  return(DT)
}
coef_to_plot_annual <- function(df, label){
  DT <- as.data.frame(df)%>%
    rownames_to_column(var = "Var") %>% 
    rename( se = `Std. Error`) %>%
    mutate( Var = str_remove(Var, ":MaxTim")) %>%
    mutate( var_n = as.numeric(str_remove(Var, "time_to_treatment::")),
            model = label,
            var_n = replace(var_n, var_n==-6, -1),
            se    = replace(se, var_n==-1, 0),
            Estimate   = replace(Estimate, var_n==-1, 0),
            ) %>%
  return(DT)
}
EV_spec <- function(df, limits, labels, colors, line_alpha){
  p <- ggplot(df, aes(x = var_n , y = Estimate, fill = as.factor(model), color = as.factor(model))) + 
      geom_hline(yintercept = 0, linetype = "dashed", size = 0.5, color = "black") +
      geom_vline(xintercept = -1, linetype = "dashed", size = 0.5, color = "black") +
      geom_errorbar(aes(ymin=Estimate-qnorm(.975)*se, ymax=Estimate+qnorm(.975)*se, color ="95% Conf."), color="grey66", width=0, size = 3.0, position = position_dodge(width = 0.7), show.legend = FALSE ) +
      geom_errorbar(aes(ymin=Estimate-qnorm(.95)*se, ymax=Estimate+qnorm(.95)*se),  width=0, size = 3.0, position = position_dodge(width = 0.7)) +
      geom_point( fill = "white", shape = 21, position = position_dodge(width = 0.7), size = 3 )  +
    geom_line(position = position_dodge(width = 0.7), alpha=line_alpha) +
    scale_color_manual(breaks = limits,
                        values = colors,
                #        limit = limits,
                        labels = labels) + 
    theme(legend.position = "none", 
          legend.title = element_blank(),
          axis.title.x = element_blank(),
          axis.title.y = element_text(size =14),
          axis.text.x = element_text(size = 14),
          axis.text.y = element_text(size = 14),
          panel.grid.major.x = element_blank(),
          panel.background = element_blank(),
          legend.background = element_blank(),
          legend.key = element_blank()
          ) + expand_limits( y = c(-1,1))
  return(p)
}


## ----MCA_Semester, include=FALSE, warning=FALSE, echo=FALSE---------------------
  setFixest_fml(..FE4 = ~ MCA + CT + AT + MT + seq)
  TIM =  
  as.data.table(read.fst("TIM_MCAS_0713.fst"))
  TIM[T== 1, effyear := seq]
  TIM[T== 0, effyear := 999999]
  TIM[, `:=`(meffyear = min(effyear),
            MaxTim = max(T),
            N = .N),
            by = list(munic, custom, country)]
  TIM[meffyear == 999999, meffyear := 0]
  TIM[, time_to_treatment := seq - meffyear ]
  TIM[MaxTim == 0, time_to_treatment := - 1000]

  ### Estimating Basic EV with Full FEs 
  EV_tim = feols(X~ i(time_to_treatment, MaxTim,
                          ref = c(-1, -1000))
                          | ..FE4, TIM, fixef.rm = 'singleton',
                          cluster = ~MCA)
  # Sun and Abraham - Estimation - Year specific Cohort
  TIM[T== 1, effyear2 := year]
  TIM[T== 0, effyear2 := 3000]
  TIM[, `:=`(meffyear2 = min(effyear2)),
            by = list(munic, custom, country)]
  TIM[MaxTim == 0, meffyear2:= 3000]

  EV_tim_SA = feols(X~ sunab(meffyear2,time_to_treatment, ref.c = 3000 )
                            | ..FE4, TIM, fixef.rm = 'singleton',
                            cluster = ~MCA)


## ----MCA_Figures----------------------------------------------------------------
# Database Coefficients
  coeftable <- coef_to_plot(EV_tim$coeftable, "1")
  coeftable_SA <- coef_to_plot(EV_tim_SA$coeftable, "1")
# Replies Figure
  ggthemr("flat")
# Paper Figure - 4
  limits <- c("1")
  labels <- c("MDC + ST + DT + CT + MT")
  colors <- c("grey30")
  ggthemr("flat")
  q <- EV_spec(coeftable%>%filter(model==1), limits, labels, colors, 0.0)
  name <- paste(FigFolder, "EV_MDC_Semester_Paper_B2.png", sep = "")
  ggsave(name, q, dpi=600, height = 5.78, width = 10.5)
  # Sun & Abraham Approach - Figure A.3
  r <- EV_spec(coeftable_SA%>%filter(model==1), limits, labels, colors, 0.0)
  name <- paste(FigFolder, "EV_MDC_Semester_Paper_SA.png", sep = "")
  ggsave(name, r, dpi=600, height = 5.78, width = 10.5)


## ----MCA_Annual, include=FALSE, warning=FALSE, echo=FALSE-----------------------
  setFixest_fml(..FE4 = ~ MCA + CT + AT + MT )
  TIM =  
  as.data.table(read.fst("TIM_MCAS_0713.fst"))
  TIM[T== 1, effyear := year]
  TIM[T== 0, effyear := 999999]
  TIM[, `:=`(meffyear = min(effyear),
            MaxTim = max(T),
            N = .N),
            by = list(munic, custom, country)]
  TIM[meffyear == 999999, meffyear := 0]
  TIM[, time_to_treatment := year - meffyear ]
  TIM[MaxTim == 0, time_to_treatment := - 1000]

  ### Estimating Basic EV with Full FEs and Standard FEs.
  EV_tim = feols(X ~ i(time_to_treatment, MaxTim,
                          ref = c(-1, -1000))
                          | ..FE4, TIM,
                          cluster = ~MCA)
  # Sun and Abraham - Estimation - Year specific Cohort
  TIM[T == 1, effyear2 := year]
  TIM[T == 0, effyear2 := 3000]
  TIM[, `:=`(meffyear2 = min(effyear2)),
            by = list(munic, custom, country)]
  TIM[MaxTim == 0, meffyear2:= 3000]
  EV_tim_SA = feols(X ~ sunab(meffyear2,time_to_treatment, ref.c = 3000 )
                            | ..FE4, TIM,
                            cluster = ~MCA)


## ----MCA_Annual_Figures---------------------------------------------------------
# Database Coefficients
  coeftable <- coef_to_plot_annual(EV_tim$coeftable, "1")
  # Paper Figure
  limits <- c("1")
  labels <- c("MDC + DT + CT + MT")
  colors <- c("grey30")
  ggthemr("flat")
  # Figure A-2
  q <- EV_spec(coeftable%>%filter(model==1), limits, labels, colors, 0.0)
  name <- paste(FigFolder, "EV_MDC_Annual_Paper_B2.png", sep = "")
  ggsave(name, q, dpi=600, height = 5.78, width = 10.5)
  
  # Figure A-4
  coeftable_SA <- coef_to_plot_annual(EV_tim_SA$coeftable, "1")
  s <- EV_spec(coeftable_SA, limits, labels, colors, 0.0)
  name <- paste(FigFolder, "EV_MDC_Annual_Paper_SA.png", sep = "")
  ggsave(name, s, dpi=600, height = 5.78, width = 10.5)

